function [dxdt,s]=rhs(t,states,p)
    s=struct();
    dxdt=zeros(size(states));
    x=states(1);
    chi=states(2);
    thetahat=states(3);
    N=@(chi)(exp(chi^2)*cos(pi/2*chi));

    y0=p.y0(t);
    y0d=p.y0d(t);
    phi=p.phi(x);

    e=x-y0;
    Phi=phi^2/sqrt(e^2*phi^2+p.epsilon(t)^2)+y0d^2/sqrt(e^2*y0d^2+p.epsilon(t)^2);


    w=(p.c+thetahat*Phi)*e;
    u=N(chi)*w;
    dchi=p.gammax*e*w;
    dthetahat=p.sigma*Phi*e^2;

    s.u=u;
    s.N=N(chi);
    dxdt(1)=p.f(p.theta,x)+p.g*u;
    dxdt(2)=dchi;
    dxdt(3)=dthetahat;
end